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ABSTRACT 

Supermassive black holes (SMBHs) found in the centers of many galaxies are understood 
to play a fundamental, active role in the cosmological structure formation process. In 
hierarchical formation scenarios, SMBHs are expected to form binaries following the 
merger of their host galaxies. If these binaries do not coalesce before the merger with a 
third galaxy, the formation of a black hole triple system is possible. Numerical simula- 
tions of the dynamics of triples within galaxy cores exhibit phases of very high eccentric- 
ity (as high as e ~ 0.99). During these phases, intense bursts of gravitational radiation 
can be emitted at orbital periapsis, which produces a gravitational wave signal at fre- 
quencies substantially higher than the orbital frequency. The likelihood of detection of 
these bursts with pulsar timing and the Laser Interferometer Space Antenna (LISA) is 
estimated using several population models of SMBHs with masses > 10 7 M Q . Assuming 
10% or more of binaries are in triple systems, we find that up to a few dozen of these 
bursts will produce residuals > 1 ns, within the sensitivity range of forthcoming pulsar 
timing arrays (PTAs). However, most of such bursts will be washed out in the under- 
lying confusion noise produced by all the other 'standard' SMBH binaries emitting in 
the same frequency window. A detailed data analysis study would be required to assess 
resolvability of such sources. Implementing a basic resolvability criterion, we find that 
the chance of catching a resolvable burst at a one nanosecond precision level is 2 — 50%, 
depending on the adopted SMBH evolution model. On the other hand, the probabil- 
ity of detecting bursts produced by massive binaries (masses > 1O 7 M ) with LISA is 
negligible. 

Key words: black hole dynamics gravitational waves cosmology: theory pulsars: 
general 



1 INTRODUCTION 

It is well established that most galaxies host supermas- 
sive black holes (SMBHs) in their centers (|Richstone et al.l 
Il99&t ). In the past decade, compelling evidence of the 
correlation between the mass of the central SMBH and 
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the bulg e velocity dispersion and l u minosity has been col- 
lected (|Ferrarese fc MerrittJ I2OO0I; iGebhardt. et all |2000| ; 
iMerritt fcFerraresd l200ll : iTremaine et al.ll2002l h indicating 
a coevolutionary scenario for SMBHs and their hosts. On a 
cosmological scale, galaxy formation and evolution can be 
understood by semi-analytic modeling, where properties of 
the baryonic matter are followed in the evolving dark matter 
halos obtained from large-scale models of hierarchical gravita- 
tional structure formation. A simple model of galaxy and cen- 
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tral SMBH evolution in which every merger of galaxies leads 
quickly to coalescence of their central black holes can quan- 
t it at ivel y reproduce both the SMBH mass- bulge luminosity 
relation (Kauff mann fc Haehneltll2000l) and the SMBH ma ss- 
velocity dispersion relation |Haehnelt fc Kau ffmann 2000). 

In this general picture, if both of the galaxies involved 
in a merger host a SMBH, then the formation of a SMBH 
binary is an inevitable stage of the merging process. Fol- 
lowing the merger, the two black holes sink to the cen- 
ter of the merger remnant becau se of dynamical friction 
ijBegelman. Blandford fc~R ees 1980|). When the mass (either 
in gas or stars) enclosed in their orbit is of the order of their 
own mass, they start to feel the gravitational pull of each 
other, forming a bound binary. The subsequent binary evolu- 
tion is, however, still unclear. In order to coalesce, the binary 
must shed its binding energy and angular momentum; a dy- 
namical process known in literature as 'hardening'. A crucial 
point in assessing the fate of the binary is the efficiency with 
which it transfers energy and angular momentum to the sur- 
rounding gas and stars. 

The case of SMBH binaries in stellar environments has 
received a lot of attention in the last decade. The system is 
usually modeled as a massive binary embedded in a stellar 
background with a given phase space distribution. The re- 
gion of phase space containing stars that can interact with 
the S MBH binary in one orbital period is known as the loss 
cone ( Frank & Rccs 197t j; lAmaro-Seoane fc Spurzeml l200ll ; 
iMilosavlievic fc Mer ritt 2003). As the binary evolves, it ejects 
stars on intersecting orbits via the so called 'slingshot mecha- 
nism', causing a progressive emptying of the loss cone, which 
ultimately increases the hardening time scale. Without an ef- 
ficient physical mechanism for repopulating the loss cone, the 
binary will never proceed to small separations where coales- 
cence induced by gravitational radiation takes place within 
a Hubble time. This is known as the st alling or 'last parsec' 
problem (|Milosavlievic fc Merrittll200lh . 

In the last decade, several solutions to the stalling 
issue have been proposed. Axisymmetric or triaxial stel- 
lar distrib utions ma y significantly shorten t he coalescence 
times cale (|Yu1 |2002| ; iMerritt fc Poonl 12004 IBerczik et~ai] 
120061 ). This is bacause the presence of deviations from 
spherical symmetry can produce "boxy" orbits, as seen by 
IBerczik. et al\ (|2006h . These orbits produce centrophilic stel- 
lar orbits and, therefore, repl enish the loss-cone. However, 
more recent calculations by lAmaro-Seoane fc Santamarial 
HH) of the outcome of the merger of two clusters ini - 
tially in parabolic orbits l|Amaro- Seoane fc Freitad [2006) 
have not been able to reproduce the rotation necessary to 
create the unstable bar structure. Other studies have in- 
voked eccentricities of the binary to refill the loss cone, 
since this effect could alter the cross section for super- 
elastic scatterings (thus altering the state of the loss cone) 
and shorten the gap to the onset of gravitational radia- 
tion effects (e.g.: iHemsendorf. Sigurdsson. fc Spurzeml [20021; 
Aarseth 2003a: Berc zik. et aZ.l 2006: Am aro-Seoane fc Freitaa 
20061 : lAmaro-Seoane. Miller fc Freitag||2009l ). The presence of 



ble quick coalescence (jEscala et alj|2005l : iDotti et alj|2006af ). 
However, current simulations do not have the resolution to 
follow the binary fate down to the gravitational wave (GW) 
emission regime, and robust conclusions about its late inspiral 
and coalescence can not be drawn. In any case, very massive 
low redshift systems, which are the major focus of our study, 
are more likely to reside in massive gas poor galaxies and 
their dynamics is probably dominated by stellar interactions. 

When scaled to very massive binaries (masses > 1O 8 M0), 
the inferred coalescence timescales in a stellar dominated en- 
vironment are of the order of few Gyrs, indicating that SMBH 
binaries may be relatively long living systems. If the typical 
timescale between two subsequent mergers is comparable the 
SMBH binary lifetime, then a third black hole may reach 
the nucleus when the binary is still in place, and the forma- 
tion of SMBH triplets might be a common step in the galaxy 
formation process. Recent studies of galaxy pairs lead to the 
conclusion that 30 — 70% of present day mass ive galaxies have 
undergone a m ajor merger since redshift one l|Bell et al.ll2006l : 
iLin et al]|2008h . where 'major' means with baryonic mass ra- 
tio of the two components larger than 1/3 or 1/4 (depend- 
ing on the study), which is a quite conservative threshold. 
This means that, on average, all massive galaxies have expe- 
rienced a merger event in the last ten billion years. Assuming 
uncorrelated events, and a typical binary lifetime of one bil- 
lion years, then 10% of SMBH binaries may form a triplet. 
With increasing redshift (and decreasing masses), dynamical 
timescales become shorter and shorter, implying that triplets 
may have been more common in the high redshift Universe. 

In this paper we focus on SMBH triplets, studying 
their dynamical evolution, GW emission, and detectabil- 
ity. Employing sophisticated three body scattering experi- 
ments calibrated on direct-summation Nbody simulations, 
we study the dynamical evolution of the system, finding 
surprisingly high eccentricities of the inner SMBH binary 
(up to e > 0.99). Even though the triple interaction would 
aossibly lead to an ejection of one or even all SMBHs 



possib ly lead to an ejection ol one or even all blvllaris 
(jValtonen. et al\ 1994), most of t he sy stems are long liv- 
ing (~ 10 9 yrs, iHoffman fc Loebl l|2007l) ). and final coales- 
cence is more common than ejection , confirming analytical 
results bv lMakino fc Ebisuzak i ( 1996). W e mod el at the lead- 
ing quadrupole order (|Peters fc Mathewsj|l963T ) the bursts of 
gravitational radiation emitted in the highly eccentric phase, 
assessing detectability with future GW experiments. Adopt- 
ing cosmologically and astrophysically motivated models for 
SMBH formation and evolution, we estimate reliable event 
rates. 

In order to cover the low frequencies generated by 
the expected c osmological population of coalescing SMBH 



binaries (e.g.. | Wvithe fc Loeb I 120031: iSesana et al. I 1200 
2005; Scsa na. Volonteri fc Haardt 1 120071 ') or plunges of com- 
pact objects such as stellar black holes on to supermas- 
sive ones (see e.g. lAmaro-Seoane et al.l 120071 . for a review 
and references there in) , the space-born observatory LISA 
(|Bender et al. I [l~998) has been planned to be covering the 
range of frequencies of ~ 10~ 4 — 10 _1 Hz. Moving to 
even lo wer frequencies, the Pa rkes Pulsar Timing Array 
(PPTA, lManchester1l2006l. 120081). th e Euro pean Pulsar Tim- 
ing Array (EPTA, Ijanssen et all 120081 ) and the North 
American Na nohertz Observato ry for Gravitational Waves 
(NANOGrav, Ijenet et alj|2009h are already collecting data 
and improving their sensitivity in the frequency range of 



massive perturbers may also help re plenishing the loss cone, 
boosting the binary hardening rate (jPerets et al.ll2007l ). On 
the other hand, in smooth particle hydrodynamics simula- 
tions of SMBH binaries in gas-rich environments, efficient 
hardening induced by the tidal interaction between the binary 
and the gas medium has been observed, indicating a possi- 
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~ 10 8 — 10 6 Hz, and in the next deca de the planned Square 
Kilometer Array (SKA, lLazio 1 120091 ) will provide a major 
leap in sensitivity. 

Throughout this paper we consider only very massive 
systems, with total mass ~ 1O 8 M0. Our goal is to in- 
vestigate if the high frequency nature of eccentric bursts 
can provide information about systems which would oth- 
erwise emit outside the frequency windows of the planned 
GW experiments quoted above, by shifting wide (sepa- 
ration > 0.1 pc) SMBH binaries into the PTA window 
or by boosting relatively massive (masses > 1O 7 M0) sys- 
tems into the LISA domain. We note that the bursts an- 
alyzed here are different from the 'bursts with memory', 
which arise during t he actual coalescence of SMBH b i naries 
and are discussed in IPshirkov. Baskaran fc Postnov I {2009) 
and Ivan Haasteren fe Levin I l|2009h . 

The structure of the paper is as follows. In Section [2] we 
describe our comprehensive study of the dynamics of triple 
systems and investigate the eccentricity evolution of the inner 
binary by using direct-summation N— body techniques and a 
statistical 3-body sample calibrated on the TV— body results. 
In Section O we model the GW signal produced by eccentric 
bursts and we introduce observable quantities for PTAs and 
LISA. In Section 4 we construct detailed populations of emit- 
ting SMBH binaries and triplets, and we discuss our results 
in terms of signal observability and detection rates in Section 
5. Lastly, we briefly summarize our results in Section 6. 



2 DYNAMICS OF TRIPLE SYSTEMS 

In modeling the dynamics of black hole triple systems within 
the centres of galaxy merger remnants, direct A r -body inte- 
grations provide the most accuracy but are the most com- 
putationally expensive. We performed eight direct TV-body 
calculations and used these to test the validity of an ap- 
proximation scheme involving three-body SMBH dynamics 
embedded in a smoothed galactic potential with dynamical 
friction and gravitational radiation modeled by drag forces. 
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Figure 1. Initial conditions for the different direct iV-body sim- 
ulations in the Hp, V plane. For each simulation we choose the 
separation between two SMBHs to be substantially smaller than 
the distance to the third SMBH. The initial parameters are se- 
lected in such a way that they are random but with an initial 
velocity smaller than the escape velocity V e sc and with a radius 
r < Rp. We also show the circular velocity in the figure, V^j rc . 
Initially The set of 3 SMBHs of simulation A is represented with 
red solid bullets; for simulation B with green solid triangles; for 
simulation C with cyan open triangles; for simulation D with blue 
open squares; for simulation E with pink open circles; for simula- 
tion H with solid orange squares and for simulations F and G with 
black crosses. In the case of simulations A, B, C, D, E, F and G 
the three SMBHs are set initially in a planar configuration. In the 
case of simulation H they have a z component different from zero 
in both the coordinates and velocities. In the cases of simulations 
A, B, C, E and H we slightly modified the positions of the symbols 
in order to avoid an overlap 



2.1 Direct TV— body calculations 

The direct-summation Nbody method we employed for all 
the calculations includes the KS regularisation. Thus, when 
two particles are tightly bound to each other or the separa- 
tion between them becomes very small during a hyperbolic 
encounter, the system becomes a candidate to be regularised 
in order to avoid problematic al small individual time steps 
(|Kustaanheimo fc Stiefel 19651) . This procedure was later ex- 
ported to systems involving more than two particles. In 
particular, the KS regularisation has been adapted to iso- 
lated and perturbed 3- and 4-body systems — the so-called 
triple (unperturbed 3-body subsystems), quad (unperturbed 
4-body subsystems) and the chain regularisation. The latter 
is invoked in our simulations whenever a regularised pair has 
a close encount er with another single star or another pair 
(|Aarsethll2003bl ). 

The basis of direct Nb ody codes rel ies on an improved 
Hermite integrator scheme (|Aarsethlll999l ) for which we need 
not only the accelerations but also their time derivative. The 
computational effort translates into accuracy so that we can 
reliably keep track of the orbital evolution of every particle 



in our system. In order to make a highly accurate estimate 
of the eccentricity evolution of the SMBH system, we do not 
employ a softening to the gravitational force (i.e. substituting 
the 1/r 2 factor with l/(r 2 + e 2 ), where r is the separation 
and e the softening parameter) that weakens the interaction 
at small separations. 

The initial conditions for the set of three SMBHs, used 
to conduct an exploration of the initial parameter space are 
shown in Table 12.11 F or the stellar system, we use a Plum- 
mer model l|Plummer|[l9Tll ). which is an n = 5 polytrope 
with a compact core and an extended outer envelope. In this 
model the density is approximately constant in the centre and 
drops to zero in the outskirts, <j> — —GM+/ \/ r 2 + R P , with 

the total stellar mass. This defines the Plummer radius 
Rp . We depict the initial conditions in Figure [T] relative to 
the circular and escape velocity of the Plummer potential. We 
present results from eight direct numerical simulations, one 
using 512,000 stars using the special-purpose GRAPE6 sys- 
tem and the remaining simulations using Beowulf PC clusters 
and the AEI mini-PCI GRAPE cluster Tuffstein. 
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Model 


A 


B 


c 








G 


H 


A/"* 


64,000 


64,000 


64,000 


64,000 


64,000 


64,000 


512,000 


256,000 


R/Rp 


0.40012 
0.40012 
0.40012 


0.80006 
0.20303 
0.20303 


0.40012 
0.64031 
0.40012 


0.80006 
0.53956 
0.20303 


0.20025 
0.20303 
0.20303 


0.60017 
0.64031 
0.40012 


0.60017 
0.64031 
0.40012 


0.49497 
0.20278 
0.20278 


V / Vesc 


U.U / 44 / D 

0.07476 
0.07476 


U.U / Q 1 

0.07476 
0.07476 


U.U / *± 1 

0.07476 
0.07476 


u. / uuuu 
0.09345 
0.09345 


0.07476 
0.07476 


fl fiAQQI 

u.o^tyy i 
0.07476 
0.07476 


u.o^tyyi 
0.07476 
0.07476 


U. tut 11 

0.09345 
0.09345 


v/v ciTC 


0.20886 
0.20886 
0.20886 


0.13543 
0.37956 
0.37956 


0.20886 
0.15108 
0.20886 


1.26802 
0.20886 
0.47442 


3.84516 
0.37956 
0.37956 


1.36391 
0.15108 
0.20886 


1.36391 
0.15108 
0.20886 


1.68376 
0.47496 
0.47496 



Table 1. Initial conditions for the set of three SMBHs in each of the eight direct TV-body simulations. A/* is the total number of stars 
employed in the simulation, V their velocity, R their position, Rp the position of the SMBHs in terms of the Plummer radius, V e sc their 
escape velocity and V c i rc is their circular velocity. The mass of the SMBHs in N— body units is 1 • 10 -2 , they are equal-mass and the mass 
of a star 1.15 ■ 10" 5 



2.2 Three-body improved statistics 

While direct iV-body simulations yield a very accurate result, 
they should be seen as a way to calibrate and test faster, 
more approximate simulations which can exhaustively cover 
the parameter space and provide good statistics. We note 
that the SMBHs in the TV— body simulations are equal-mass 
and (with the exception of simulation H), all of the systems 
studied with this method are coplanar. This was done because 
setting all SMBHs on a single plane accelerates the dynamics, 
shortening the integration time. 

In general, one wants to explore the whole parameter 
space, including non coplanar systems with different SMBH 
masses. For this purpose, we performed an ensemble of 1000 
three-body experiments, with the three Euler angles of the 
outer orbit sampled uniformly and a distribution of mass 
ratios motivated by Extended Press-Schechter theory (with 
typical mass ratios mi : m2 : m3 around 3.5:1). In each exper- 
iment we computed the Newtonian orbits of three SMBHs 
embedded in a smooth galactic potential and added drag 
forces to account for gravitational radiation and dynamical 
friction. We also included coalescence conditions when ei- 
ther of the two SMBHs pass within three Schwarzschild radii 
of each other, or the gravitational radiation timescale be- 
comes short relative to the orbital period of the binary. Close 
triple encounters were treated using a ICS-regularised few- 
body code provided by Sverre Aarseth (jMikkola fc Aarsethl 
Il990l . ll993h . while the two-body motion in between close en- 
counters was followed with a simple 4 th-order Runge-Kutta 
integrator. See lHoffman fc Loebl (|2007T l for further details on 
the code. The initial condi t ions a re those for the canonical 
ICs as in iHoffman fc Loebl (|2007T j . We performed each run 
twice — once with gravitational radiation drag and the coa- 
lescence conditions, and once without. 

The 3-body experiments are divided into two computa- 
tional regimes based upon a dimensionless parameter, a, that 
measures the relative tidal perturbation to the inner binary 
by the interloper at apoapsis: 

„ -Rapo -^single , . 

a = 2 Tf TP ' W 

^"bin, smaller -'-single 

where R apo is the apoapsis separation of the two inner binary 
members, M s i n gic is mass of the single SMBH, Mbin, smaller is 



mass of the smaller inner binary member and Dsmgie is the 
distance of the interloper (single SMBH) from the inner bi- 
nary centre-of-mass. In the limit a — » 0, we know that the 
period of the inner binary is perfectly Keplerian (plus gravita- 
tional radiation) , since the perturbation to the force from the 
third body is negligible, and thus we can do orbit-averaged 
integration instead of precisely integrating the trajectories of 
all three bodies. The two regimes are defined as follows: 

(i) The first regime corresponds to when a 3-body inter- 
action is taking place (denned by a > 10 -5 ), an extremely 
conservative criterion for when we need to do the full three- 
body integration. 

(ii) The second regime corresponds to when the single 
SMBH and binary are wandering separately through the 
galaxy (a < 10 -5 ), often on the order of a Hubble time. 

In regime (i) the 3-SMBH dynamics is integrated using Sverre 
Aarseth's high-precision, regularised CHAIN code. Gravita- 
tional radiation and stellar-dynamical friction are treated as 
perturbing, velocity-dependent forces on the three separate 
bodies. In regime (ii) the separate orbits of the single and 
binary centre-of-mass are followed using a simple 4th-order 
Runge-Kutta integrator and the evolution of the binary semi- 
major axis and eccentricity are evolved using orbit-averaged 
equations. 

In computational regime (ii) (between 3-body encoun- 
ters), the stellar int eractions are tre ated using the hard bi- 
nary prescription of iQuinlan I (|l996h . The eccentricity evo- 
lution of the inner binary under stellar interactions for near 
equal-mass hard binaries is quite weak so it is neglected en- 
tirely and only the binary semi-major axis is evolved under 
stellar interactions. The ec centricity is ev olved under gravita- 
tional radiation as given bv |Petersl(|l964l ). Dynamical friction 
tends to increase the eccentricity of a binary in the supersonic 
regime (where the orbital speed exceeds the stellar veloc- 
ity dispersion), and to circularize it in the subsonic regime. 
Since the triple SMBH system starts out supersonic in these 
3-body experiments, this effect produces a slight increase in 
the outer binary eccentricity during the initial inspiral of the 
third SMBH. Although we neglect the eccentricity evolution 
of the inner binary during this phase, we find that its eccen- 
tricity is thermalized by the first resonant 3-body encounter. 
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Figure 2. Cumulative fraction of time for the set of 1000 three-body simulations. The red, solid line corresponds to all simulations with 
the 2.5 post-Newtonian correction term; the blue, solid line corresponds to the same simulations but without it (i.e. purely Newtonian); 
the dot-dashed red curve is like the first case but taking into account only systems in which the third SMBH had an inclination below 
the critical Kozai angle of 39°, to compare with the direct Af-body simulations of Fig.Q; the dot-dashed blue curve is the same, but for 
the Newtonian cases and the dot-dashed black curve corresponds to the thermal distribution, since the direct Af-body simulations do not 
have relativistic correction terms. The top left panel shows the eccentricity computed from the instantaneous positions and velocities of 
the two binary members, and the top right panel shows the ratio of the actual to the thermal distribution. For the Newtonian runs, the 
distribution is within a factor of 2 of thermal down to 1 — e = 0.01 and within a factor of 3 of thermal down to ~ 1 — e = 0.001. The 
bottom two panels show the pericenter distances of the binary in pc and in units of the sum of the Schwarzschild radii of the two binary 
members. No coalescence was allowed during close encounters for the Newtonian runs. We can see that the distributions converge, with 
respect to the statistics, from the fact that they are substantially the same as the lower-number experiments and that the Newtonian and 
gravitational radiation distributions match at low eccentricities 



The timescale of the chaotic 3-body interactions (~ 10 5 yr) is 
much shorter than the stellar-dynamical timescale ( 10 7 — 10 10 
yr), so any effect of stellar interactions during these encoun- 
ters is completely negligible. Thus, stellar interactions play 
only two roles in our 3-body simulations; 

a They bring the third BH in to interact with the inner 
binary from an initial marginally stable hierarchical triple 
configuration; 

b During phase (ii), stellar-dynamical interactions gradu- 
ally decrease the binary semi-major axis, so that it enters the 
next three-body encounter harder than it left the last one if 
the time between encounters is long ( > 10 7 yrs). The binary 



may even coalesce between encounters due to this gradual 
shrinking. 

Consequently, the distribution in eccentricities is best 
estimated using the fraction of time that binaries spend at 
a given eccentricity while in computational regime (i) since 
there is minimal evolution of the eccentricity while in regime 
(ii). We note that the overwhelming majority of the time 
spent in regime (i) is still spent with a small enough that the 
system can be though of as a separate inner and outer binary, 
and so the instantaneous inner binary semi-major axis and 
eccentricity are well-defined. 
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Figure 3. Time spent at different locations in the t gI - (1-e) plane in the Newtonian simulations (blue lines) of Figure 2. P = t/t c i ose 
is the probability of finding the binary in a given bin at a randomly chosen time, where t c i ose is the total simulation time spent in close 
3-body encounters. The simulation being presented here is the same as the Newtonian simulation (blue lines) in Figure 2. The purpose 
of this figure is to understand whether gravitational radiation is the main reason for the fall-off in the eccentricity distribution, since the 
highest eccentricity systems have short coalescence times and quickly disappear. f. gr is in years in the upper panel and in orbital periods 
in the lower panel (the typical resonant encounter takes on order 10 s yrs). Note that 1-e 
orbital period 



2.3 Distribution of Eccentricities 

Figured shows the fraction of the time that the binary (clos- 
est SMBH pair) spends above a given eccentricity during close 
encounters, averaged over 1000 three-body experiments. The 
red solid curves are the results of the standard runs including 
gravitational radiation drag and coalescence conditions while 
the blue solid curves are the "Newtonian" case with these 
effects neglected. The dashed red and blue curves are aver- 
aged over only those experiments where the initial BH con- 
figuration has inclination < 39°, the critical angle for Kozai 
oscillations, for comparison with the direct iV-body simula- 
tions that use coplanar initial conditions. The black (dot- 
dashed) line shows the thermal distribution of eccentricities 
for reference purposes. Three-body interactions result in a 
thermal distribution of eccentricities, truncated at very high 
eccentricities by coalescence in collisions when gravitational 
radiation is included. The similarity of the dashed and solid 
curves shows that once the initial secular evolution is over, 



the system quickly thermalizes and little memory of the ini- 
tial configuration is maintained. Figure 05 shows the ratio 
of the thermal to the actual distribution as a function of ec- 
centricity. The first close encounter in each experiment has 
been excluded from these plots, since it begins from a stable 
hierarchical triple configuration and includes a long period of 
secular evolution, whereas chaotic three-body encounters are 
the focus of this work. 

The runs with and without gravitational radiation 
closely follow each other and the thermal distribution up to 
eccentricities e ~ 0.99. At higher eccentricities, the Newto- 
nian distribution remains within a factor of 2 to 3 of the 
thermal distribution, but the gravitational radiation curve 
falls off sharply, since these high-eccentricity systems coalesce 
quickly through emission of gravitational waves. To verify this 
interpretation of Figure [5] we plot the time spent at different 
locations in the t gr -(l — e) plane in Figure [3] where the grav- 
itational radiation timescale is computed from the instanta- 
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neous (a, e) of the binary. The eccentricity where the red and 
blue curves diverge in Figure [2)3 (1 — e ~ 0.01) is the value 
where the typical t gr falls to just a few orbital periods, so that 
the binary can coalesce quickly by gravitational radiation be- 
fore the third body scatters it on to a lower-eccentricity or- 
bit. Figures 0>d show the time spent by the binary at various 
pericenter separations, in pc and in Schwarzschild radii. Note 
that the gravitational radiation curve diverges sharply from 
the Newtonian one when the pericenter separation reaches 
~100 Schwarzschild radii. We show the fraction of time spent 
at different eccentricities and the fraction of time spent at dif- 
ferent pericenter separations for the JV— body simulations in 
Figure [4] The qualitative features are retained. 



3 GRAVITATIONAL WAVES: ANALYSIS OF 
THE SIGNAL 

In this section we make use of the leading Newtonian order 
derivation of the GW radiation from eccentric binaries, as 
described in iPeters fc Mathews! (|l963l ). We also use of ge- 
ometric units, with G = c = 1. Consider a system with 
masses M2 < Mi orbiting with an orbital rest frame fre- 
quency f r = uj/2tt, and with eccentricity e; at the quadrupolc 
leading order, the luminosity E emitted by the system aver- 
aged over one complete orbit is: 



32 



E = ^-M w/ "{2nf r ) w/ "F(e) = E c F(e) (2) 



where 



f(e) = X>(n,e) = — ~* 



1 + 23 e 2 , 37 4 
T 24 c T 96 c 



(1 - e 2 ) 7 /2 



(3) 



and M = Ml' 6 Ml 15 /{Mi + Af 2 ) 1/5 is the chirp mass of the 
system. E c , defined by the right-hand side in equation ([2]), is 
the luminosity emitted by a circular binary orbiting at the 
same frequency f r . The binary radiates GWs in the whole 
spectrum of harmonics / r , n = nf r (n — 1,2,...), and the 
relative power radiated in each single harmonic is described 
by the function g(n, e), defined as: 



9(n,e) = 2L (b* + (1 - e 2 ) A 2 n + JLj n (ne) 2 ) , (4) 

where J„(x) are the Bessel functions and A n and B n are also 
defined in terms of the J n as: 



B n = Jn-iine) — 2eJ„-i(ne) H — J„(ne) 

n 

+2eJ„+i(ne) - J n+2 (ne) 
A n = J„- 2 (ne) — 2J„(ne) + J„ +2 (ne). 



(5) 
(0) 



The total luminosity of the source can then be written, using 
equations ((2} and as the sum of the component radiated 
at each single harmonic: 



(7) 



n=0 n=0 



Given a general GW characterised by the two polarised 
component waves h+ and h x , the rms amplitude of the wave 

is defined as h = (h 2 ^ + /i 2 }, where ( ) denotes the average 
over directions and over time. The flux radiated in the GW 



field is related t o the derivativ es of its amplitude components 
by the relation (|Thorndll987T ) 



dE 
oltdA 



(Si 



The sinusoidal nature of the waves implies {h\ + ft 2 } = 
47r 2 / 2 {h\-\-h 2 x ) . So that integrating equation ([8]) over a spher- 
ical surface of radius &l (the luminosity distance from the 
source) centered at the source and averaging over an orbital 
period, directly relates E to the wave rms amplitude. We can 
then infer that the rms amplitude and the energy radiat ed in 
the n-th harmonic are related as (|Finn fc Thornd feoOO) 



h n = 



Tld f r . 



= 2 



32 M 



5/3 



5 ndr 



(2nf r ) 2 ^^M, (9) 



where d = di/(l + z). In the limit of a circular orbit (i.e., 
g(n,e) = <5 n ,2 in the Kronecker-5 notation), equation re- 
turns the usual sky-polarization averaged amplitude l|Thornel 
Il987t l. 



3.1 Observed quantities 

Since we are interested in an estimate of the detectability of 
extremely eccentric binaries (induced by triple interactions) 
by means of pulsar timing (and possibly LISA) observations, 
we first introduce an extension of the characteristic amplitude 
to include eccentric binaries. Eccentric binaries emit pulses 
of GWs at their periapsis passages, and the rms amplitude 
of each harmonic is given by equation ([9]|. However, h n is an 
average amplitude related to the average luminosity along the 
orbit. The actual relevant time for the burst is the periapsis 
passage timescale T p = (1 — e) 3//2 T or b (T or b is the binary 
orbital period), and if the burst is detected, almost all the 
energy radiated along the whole orbit is seen on the timescale 
T p . This means that the relevant detectable amplitude of each 
harmonic during the burst is 



^obs,n — hn y/ 7 ' fr 



(10) 



where the factor T = max(T or b, Tobs) takes into account the 
fact that, if T or b < T bs, multiple bursts are visible during the 
observation. Equation (|10|l is a crude approximation, never- 
theless it catches the basic features of the observed signal: 
this is given by the rms amplitude of each single n-th har- 
monic multiplied by the square root of the cycles completed 
by the harmonic in an orbital period, assuming that the bi- 
nary orbit is a fixed ellipse and GW emission does not change 
the orbital parameters. 

The search for GWs using pulsar timing data exploits 
the effect of gravitational radiation on the propagation of 
the radio waves from one (or more) pulsar(s). A pass- 
ing GW would imprint a characte ristic s i gnature on the 
time of arrival of radio pulses (e.g . ISazhinl Il978l : iDetweilerl 
1 19791 : iBertotti. Carr fc Reesl Il983h . pro ducing a so called 
tim ing residual. We refer the reader to Ijenet et al.l l|2004h 
and lSesana. Vecchio fc Volonteril <|2009l ), (hereinafter SVV09) 
for a detailed mathematical description of the GW induced 
residuals. The residuals are defined as integrals of the GW 
during the observation time. For a collection of harmonics, 
the residuals are given by: 



R(T) 



E 



p 2 



2(1+7) 



a;3 



(1+7) 



dt, (11) 
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Figure 4. Left panel: Fraction of time for the direct TV-body simulations at which the binary black hole has eccentricity in the range 
1 — e for the 512,000 stars simulation (thin green curve) and for the lower-resolution simulations. The thick red curve corresponds to the 
average of these lower-resolution computations and the dashed curve to the thermal distribution. Right panel: Cumulative fraction of time 
spent on a certain periapsis distance for all direct TV-body simulations following the same colour labelling. One TV— body unit of distance 
is U\n = 1.1 pc 



where a, j3, and 7 are the direction cosines of the pulsar 
relative to a Cartesian coordinate system defined with the 
z-axis along the direction of propagation of the gravitational 
wave and the x and y axes defining the + polarization. The 
harmonics of the two polarizations , h+, n and h x ,n, can be 
found in Section 3.2 of lPierroT et aZ.I(|200ll) . The rms residual 
5t S w is then formally defined as w (R(T) 2 ). 

A simple derivation of the average timing residual 5t gw 
generated by a circular binary is given by SVV09. With the 
notations adopted above, their equation (20) reads: 



and cosines of different harmonics, that drop to zero when av- 
eraged over the observation, leaving only a sum of the square 
signals produced by each single harmonic (those terms includ- 
ing cos 2 (27r/„t) and sin 2 (27r/„£)). We shall plot, in Section 4, 
R(T) for selected eccentric bursts, and we will see that <5t gw 
as defined by equations 1)130 and 1 )140 gives a good estimate 
of the amplitude of the induced residual. 

For inferring LISA detectability, given /i bs,m an estimate 
of the signal to noise ratio (SNR) in the LISA detector is 
straightforwardly computed as: 



ftgwCf) = 



15 27T/ r 



(12) 



where the observed frequency / is related to f r as / = / r /(l + 
z) (being z the redshift of the source), the factor \//T b s takes 
into account for the signal 'build-up' with the square root 
of the number of cycles, and y8/15 comes from the angle 
average of the amplitude of the signal (cf. equation (17)-(21) 
of SVV09). We can generalise this derivation to the case of 
bursts produced by eccentric binaries, relating the h bs,n of 
each harmonic to the induced residual residual at its peculiar 
frequency via: 



^gw(/n) — 



8 h 



(13) 



15 2Tvf r ,n 

The total residual can then be assumed to be of the order: 
/« \V> 

*tgw= £*&r(/n) • (14) 

\n=0 / 

The estimation given in equations (|13|l and ()14[) is justified 
because the integral in equation (|ll)l gives products of sines 



SNR 2 = 4 J2 



5/5/ 



(15) 



where Sf is the one-side noise spectral densi ty of the detector. 
We ad opted the Sf given in equation (48) of lBarack fc Cutler! 
(|2004h , based on the LISA Pre-Phase A Report. We extended 
the sensitivity down to 10 -5 Hz and we considered detec- 
tion with two independent TDI interferometers (which im- 
plies a gain of a factor of two in Sf). The SNR computed in 
this way may seem a poor approximation. However, we have 
checked the SNRs against th ose obtained follow i ng th e proce- 
dure given in Section V-B of lBarack fc Cutled l|2004h . where 
the binary is consistently evolved with orbit averaged post- 
Newtonian equations, and found agreement at a 20-30% level, 
which is acceptable since we are interested in a preliminary 
estimation of source detectability 0. 



1 The difference is mainly due to the fact that the orbital param- 
eters change during the strong GW emission burst, and this is not 
taken into account in equations dlOI I and 1151 . 
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3.2 Some heuristic considerations 

The previous derivation can be use to achieve a heuristic 
understanding of what we may expect to actually detect. Let 
us consider two binaries '1' and '2' with the same masses, and 
semimajor axes related as ai = ai(l — e) (suppose '2' is in 
circular orbit and '1' on a very eccentric orbit, i.e. 1 — e <C 1). 
Equations ([2} and (£3j) provide the luminosity averaged over an 
orbital period. The eccentric binary '1' has an orbital period 
T\ oc a 3 / 2 . But, it emits GWs in a short burst of duration 
of the order of its periapsis passage that is T p tx [ai(l — 
e)] 3 / 2 . The mean luminosity of the eccentric binary during 
the periapsis burst is then 



pTi F(e) 3/2 
^iTTT cx —5-1,1 - e) 

1 p d-^ 



(16) 



where we used the Newtonian relation to switch from / to a 
in equation ([2|, and we ignored the source redshift. According 

to equations (|8]) and ((9|, we can write hi ~ \J Ei, P / f P , where 
we make the assumption that f p is the 'dominant frequency 
of the burst', which corresponds to the 'periapsis frequency', 
f p ~ /i(l — e)~ 3 ^ 2 . A circular binary with semimajor 02 

simply emits a periodic wave with amplitude /12 ~ \/ .E2//2, 
where E2 oc l/a\. Remembering that ai — ai(l — e) and, 
consequently, f p ~ f, the /11//12 ratio reads: 



hi 
hi 



1 + 73 2 , 37 4 

'24 ~ 96 

(l + e) 7 /2 



1/2 



0(1). 



(17) 



Since PTAs detect a timing residual that is St ~ ft-//, it 
follows that Sti/8t2 ~ hi/h2- The timing residual caused by 
a burst that happens to be at the right frequency for PTA 
(~ 10~ 8 Hz), generated by a very eccentric binary with an 
orbital frequency / <C 10 -8 , is then of the same order of the 
residual caused by a circular binary emitting at / = 1CP 8 . 
The signal is, however, quite different and it is spread over a 
broad frequency band. This heuristic consideration suggests 
that PTA detection of such extreme events may be rather 
difficult, because their signal may be overwhelmed by GW 
emitted by 'conventional' binaries with shorter periods. On 
the other hand, we might expect some interesting effect for 
LISA, since this mechanism can boost the GW frequency 
by more than three order of magnitudes and signals from 
systems that would emit at much lower frequencies, may be 
shifted into the LISA domain. 



4 CONSTRUCTING THE SIGNAL FROM 
BINARY AND TRIPLET POPULATION 
MODELS 

4.1 Hierarchical models for SMBH evolution 

To draw sensible predictions about the number of expected 
detectable GW bursts, we need to model the population of 
triple systems that form during the SMBH hierarchical build 
up. We start by considering the SMBH binary, population. 
We are mainly interested here in probing massive systems 
M = Mi + A/2 > 10 7 Mq, so that we can use catalogs of 
syste ms extracted from the Millennium Run l|Springel et al.l 
2005). We employ the very same catalogs used in SVV09; 
the reader is referred to Section 2 of that paper for details, 



here we merely summarise the basics of the procedure. We 
compile catalogs of galaxy mergers from the semi-analytical 
model of iBertone. De Lucia fc Thomas! (|2007l ) applied to the 
Millennium Run. We then associate a pair of merging SMBHs 
to each merging pair of spheroids (elliptical galaxies or bulges 
of spirals) according to four different SMBH-host prescrip- 
tions (Section 2.2 of SVV09). Here we consider the three 
Tu models presented in SVV09, in which SMBHs correlate 
wit h the sphero i d ma sses according to the relation given 
by iTundo et all (|2007h . and differ from each other in the 
adopted accretion prescription: the Tu-SA model (accretion 
triggered on to the more massive black hole before the final co- 
alescence), the Tu-DA model (accretion triggered before the 
merger on to both black holes) and the Tu-NA model (ac- 
cretion triggered after the coalescence) . We also investigate 
the dependence on the adopted SMBH binary population by 
considering the La-SA and Tr-SA models (see SVV09 for de- 
tails). The catalogs of coalescing binaries obtained in this 
way are then properly weighted over the observable volume 
shell at each redshift to obtain the differential distribution 
d 3 N/dMdzdt r , i.e. the coalescence rate (the number of coa- 
lescences iV per unit proper time dt r ) in the chirp mass and 
redshift interval [M, A4 + dA4] and [z, z + dz], respectively. 



4.2 Signal from SMBH binaries and triplets 

The GW signal can be divided into two contributions — one 
from the binaries, and one from the triplets. We will refer 
to the latter as bursting sources, since we consider the GW 
bursts they emit at the periastron in their eccentric phase. In 
this study, we consider the binary population emitting in the 
PTA domain to be composed of circular systems dynamically 
dri ven by GW emission only. The G W signal is then given 
by (|Sesana. Vecchio fc Colacinoll2008T ): 



r 00 r 00 

h 2 c {f) = dz dM 
Jo Jo 



d s N 



dzdMdlnfr 



h 2 (fr 



(18) 



where h is the sky-polarization average of each single source 
(|Thornd 119871 ), and d s N / dzdMdlnfr is the instantaneous 
population of comoving systems emitting in a given logarith- 
mic frequency interval with chirp mass and redshift in the 
range [AI, AI + dM] and [z, z + dz], and is given by: 



d 3 A 



dzdMdlnfr 



d 3 A 



dt r 



where 



dt r 



dlnfr 647T 8 / 3 



dzdMdtr dlnfr 



5/3 f - 8/3 



M' a/A f; 



(19) 



(20) 



In equation (|19|) . Ft is the fraction of coalescing binaries that 
have experienced a triple interaction. This can be estimated 
simply by knowing the likelihood of forming triple systems 
because of two subsequent mergers. The galaxy merger rate 
drops dramatically at low redshift, and the typical timescale 
between two subsequent major merger could be as long as 
~ 10 10 yr. This means that massive galaxies may have ex- 
perienced , on average, ju st one major merger since z = 1 



perienced , on a verage, just o 
(see, e.g. iBell et al.l I2006I '). If 



we assume survival time of a 
binary is ~ 10 u yrs, adopting the simplifying assumptions 
of uncorrelated mergers with a Poissonian delay distribution 
with a characteristic time of 10 10 yr, the probability of hav- 
ing two subsequent mergers in a 10 yr time interval is ~ 0.1. 
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Figure 5. Two dimensional joint probability distribution V{r p ,e) 
for the inner binary in the [r p ,e] plane (where r p = a(l — e)). 
The distribution is obtained averaging over the 1000 3-body ex- 
periments described in Section l2.2l 



We will consider two different situations, choosing the frac- 
tion of SMBH binaries experiencing a triple interaction to be 
Tt = 0.1 or Ft = 0.5. 

By knowing Tt , we can write the coalescence rate of bi- 
naries that have experienced a triple interactions as Tt x 
d? N J ' dzdMdt r . From the 3-body scattering presented in Sec- 
tion 12.21 we derive the joint probability distribution for the 
inner binary of having a certain periastron and a certain ec- 
centricity, V(r p , e). This quantity is plotted in figure[5]for our 
set of the 1000 3-body realizations. This probability distribu- 
tion refers to systems with mean total mass ~ 4 x 1O 8 M0. 
To extend it to a wider range of masses, we assume a triplet 
lifetime T = 10 9 yrs independently of the masses (which are 
in a narrow range peaked around 1O 8 M0 in our case) and 
we rescale V(r p ,e) so that, in the GW dominated regime, 
elements in the (r p ,e) space having the same coalescence 
timescale T gw (r p , e), have the same probability value V(r p , e). 
Since T gw oc a 4 l\M\M2(M\ +M2)], assuming an invariant bi- 
nary mass ratio distribution in the relevant mass range (which 
is a good approximation given the narrow mass range we are 
dealing with), the y axis in figure [5] is rescaled for any given 
total mass of the binary M according to (M/4 x 1O 8 M ) 3/4 . 
We then compute the distribution of eccentric binaries emit- 
ting an observable burst as: 



N(M,z,r p ,e) 



d 3 N 



dzdMdt, 



x T t xT x V{r p ,e)x 
xmin[l,(T obs /r orb )]. 



(21) 



Where the factor min[l, (T b s /T or b)] takes into account the 
fact that if the binary period is longer than the observa- 
tion time, only a fraction Tobs/Torb of the systems is actually 
bursting during the observation. 



4.3 Practical computation of the signal 

The relevant frequency band for pulsar timing observations is 
between 1/T b s and the Nyquist frequency 1/(2 At) - where 
At is the time between two adjacent observations-, corre- 
sponding to 3 x 10 -9 Hz - 10 -7 Hz. The frequency resolution 



bin is Af = 1/T bs, an d we assume T b s = 10 yr through- 
out the paper. Every realistic frequency-domain computation 
of the signal has to take into account the frequency resolu- 
tion bin A/ of the observation. The signal is therefore eval- 
uated for discrete frequency bins Afj centered at discrete 
values of the frequency fj, where /y+i) = fj + A/. What 
we actually collect in our code is the numerical distribution 
A 3 N/AzAMAf r , where Af r = (1 + z)Af. The integral in 
equation (I18|) is then replaced as a sum over redshift and 
chirp mass, and the value of the characteristic strain at each 
discrete frequency fj is computed as 

A 3 N 



»c(/,) = ££ 



AzAMAf ri 



-frh\z,M,f r )AzAM, (22) 



where A/ rj = (1 + z)Afj is the jih frequency bin shifted ac- 
cording to the cosmological redshift of the sources. Equation 
(|22[) is simply read as the sum of the squares of the char- 
acteristic strains of all the sources emitting in the observed 
frequency bin Afj. If we produce a family of a = 1,...,K 
sources by performing a Monte Carlo sampling of the nu- 
merical distribution A 3 N/AzAA4 Af r of the emitting binary 
population, the characteristic strain is computed as 

K 

^(/ 3 ) = £ftc,a(2 J M ! / a , r ) 2 e[/ a , rj A/ i (l + 2)] (23) 
a = l 

where 9[/ a , r , A/, (1 + z)} = 1 if f a , r 6 A/j(l + z) and 
is null elsewhere. To recover equation (|22[1 . the character- 
istic amplitude of the individual source is given by: h 2 c a = 



hctfoL,r I Af r . 



h a fj/Afj = h a fjT ohs ; i.e., the sky and 



polarization averaged amplitude square, multiplied by the 
number of cycles completed in the observation time. The in- 
duced rms residual of each individual source is then given by 
equation (|12p . Note that in the limit of large K (formally, 
K — > 00), h c (fj) computed according to equation (|2"3")l is 
independent of T b s (because the increment of the contribu- 
tion of each single source according to the number of cycles 
completed during T b s is balanced by the fact that we sum 
over a frequency bin that is proportional to l/Tobs), and its 
value coincides with the one obt ained from the standard en- 
ergy based definition of h c {f) (|Sesana. Vecchio fc Colacinol 
l200ot ). On the other hand, when K is small (i.e. we sum over 
a small number of sources), fluctuations become important 
in the computation of the signal in each frequency bin. Nu- 
merical computation according to equation (|23[) allows us to 
account for signal fluctuations, which are missing in the an- 
alytical definit ion of the chara cteristic amplitude of the GW 
spectrum (e.g. IPhinnev i feoOl), but are important in the ac- 
tual computation of the observed signal. Given h c , the in- 
duced rms timing residual produced by the whole emitting 
population is simply given by h c (fi) / (2n fi). 

We generate a population of emitting binaries accord- 
ing to the numerical distribution A 3 N / AzAM Af r , and we 
sum all the h CtCl contributions in every frequency bin to ob- 
tain the characteristic strain of the signal. We then generate, 
again using a Monte-Carlo sampling, a population of emit- 
ting eccentric binaries in triple systems from the distribution 
given in equation (|21|l and we compute their GW bursts and 
the induced rms residuals according to equations (|9] 1101 1131 
I14|) . For the few systems reaching 10 -5 Hz with their higher 
harmonics, we also compute the SNR produced in the LISA 
detector usi ng equation l|15l). adopt ing the Sf given in equa- 
tion (48) of iBarack fc Cutlerl <|2004D , extended downward to 
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Figure 6. Representation of all the relevant features of a Monte- 
Carlo generated signal. The Tu-DA model with T = 0.5 is as- 
sumed. The jagged blue line is an individual Monte-Carlo realiza- 
tion of the signal. The small black triangles label the characteristic 
strain of the brightest source in each frequency bin. If the source 
is resolvable, it is also labeled with a big red triangle. The jagged 
red line is the stochastic level of the signal, i.e., once the resolvable 
sources in each frequency bin are subtracted. The magenta points 
label all the systems producing an rms residual (computed through 
equation i!2l ) larger than 0.1ns over 10 years. The cyan 'arcs' of 
dots, represent the contribution to the signal coming from eccen- 
tric binaries in triple systems (bursting sources), again assuming 
that their total rms residual is larger than 0.1ns (equation J 141 1. 
The arc-like black (red) tracks represent the spectrum of the more 
luminous (resolvable) bursting systems in the realization, and have 
the only purpose of guiding the reader eye. The dotted oblique lines 
mark different rms residual levels as a function of the frequency. 

10~ 5 Hz as described in Section f3.ll We consider five differ- 
ent SMBH binary populations presented in SSV09 (Tu-SA, 
Tu-DA, Tu-NA, La-SA, Tr-SA) with two different fractions 
of triplets Tt = 0.1, 0.5, for a grand total of 10 different mod- 
els. We run 50 (when Tt = 0.5; 100 if Tt = 0.1) independent 
Monte-Carlo realizations of each single model, which allows 
us to perform a statistical study of the properties of the burst- 
ing sources. We consider only systems with e > 0.66, because 
highly eccentric systems are those expected to burst at high 
frequencies, where the contribution of the overall circular bi- 
nary population declines. And also because high eccentricities 
result in a well defined burst shape which may be essential 
to distinguish it from periodic sources. 



5 RESULTS 

5.1 description of the signal 

All the relevant features of the signal are plotted in figure 
[5] for a realization of the Tu-DA model with Tt = 0.5. A 
Monte-Carlo generated signal is depicted as a blue jagged 
line. The magenta points represent all the binary systems 



producing a St gv > 0.1 ns; there are ~ 4000 sources in this 
particular realization. The cyan 'arcs' of dots, represent the 
contribution to the signal coming from eccentric binaries in 
triple systems (bursting sources), where contributions from 
all harmonics falling in the same frequency bin were added 
in quadrature. The black triangles correspond to the bright- 
est source in each frequency bin. And if a source is brighter 
than the sum of all the contributions coming from the other 
sources emitting in the same bin, we consider that source 
resolvable and we mark it with a superposed red triangle. 
The red jagged line is the resulting stochastic level of the 
signal, after the contribution from the resolvable sources has 
been subtracted. The arc-like black (red) tracks represent the 
more luminous (resolvable) bursting systems in the realiza- 
tion. In this particular case there were five resolvable bursts 
with rms residual 8t gw = 3.5,0.07,0.04,0.01,0.002 ns. How- 
ever, considering realistic PTA sensitivities achievable in the 
near future (~ Ins, with the SKA), only the brightest one 
would have a good chance of being detected. 

We note that we introduced the concept of resolvable 
source in the frequency domain, assuming that a source is 
resolvable if its strain is larger than the sum of the strains 
of all the other sources in that frequency bin. This definition 
is, however, only appropriate for monochromatic sources. A 
very eccentric burst, emitting a whole spectrum of harmonics, 
may not be the brightest source in any of the frequency bins, 
however, it may produce a significantly larger rms residuals 
with respect to other individual circular binaries. Moreover, 
in the time domain, the signature of these bursts is quite 
different with respect to periodic circular binaries, resulting 
in long bumps or narrow well localized bursts (see figure [10}. 
Given these caveats, we will also present results in terms of 
total number of sources, independent of their resolvability 
according to our definition. 

A sample of different realizations of the signal is col- 
lected in figure [7] for individual realizations of the three dif- 
ferent Tu models. Only the brightest sources are plotted in 
this case. Given the small number of systems involved, their 
phenomenology is quite variable. For example, the realiza- 
tion illustrated in the left-middle panel, shows three resolv- 
able bursts with 8t gv > 3ns; the one in the lower- left panel, 
does not show any individually resolvable bursts. 

5.2 Statistic of bursting sources 

To quantify the statistics of the bursting sources, we cast the 
results in terms of the cumulative number of sources as a 
function of the timing residuals: 

N(5t g „) = / — — -d(St' gvJ ) , (24) 

where the distribution dN/d(5t' gw ) is the average over the 
50 (100) Monte-Carlo realizations of each model. We com- 
pute this average both considering all the sources emitting 
over a given <5t gw threshold (obtaining the total distribution 
of bursting sources), as well as considering only resolvable 
sources as defined in the previous section (obtaining the dis- 
tribution of bursting resolvable sources). In figure [8] N(St gw ) 
for all the sources is shown. Depending on the adopted model, 
and on the fraction of triplets assumed, there are few hun- 
dred to few thousand binaries contributing to the signal at 
a level > Ins. The number of triplets over this threshold is 
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IO" 8 10- 7 io- 8 10- 7 io- 8 10- 7 

observed frequency [Hz] 

Figure 7. Sample of individual Monte— Carlo realizations of the signal generated using Tu-DA (left panels), Tu-SA (central panels) and 
Tu-NA (right panels) models (J 7 = 0.5). Line and point style as in figure[6] The two dashed lines in each panel represent the sensitivity of 
the PPTA (upper) survey and an indicative sensitivity of Ins for SKA (lower). 



between 20 and 60 assuming T = 0.5, and, not surprisingly, 
a factor of 5 lower if we assume T = 0.1. If triple interactions 
of SMBHs are common (say, T > 0.1), we may therefore ex- 
pect l-to-100 bursts from eccentric sources contributing to 
the GW signal at a residual level of > Ins. The eccentricity 
distribution of these bursts is basically fiat in the considered 
eccentricity range (0.66,1). If we consider resolvable sources 
only, the figures are not as promising. As shown in figure [9] 
a timing precision of 0.1 ns is needed to guarantee the de- 
tectability of a resolvable burst if T = 0.5. At a 1 ns level, 



we have less than one resolvable burst, we can then interpret 
the results in terms of the probability of having such bursts 
in our observable Universe. This probability ranges from 2% 
to 50% depending on the adopted model, and the eccentric- 
ity distribution of these resolvable events is biased towards 
high values, peaking around e = 0.9. La-SA and Tr-SA give 
similar results both qualitatively and quantitatively, there- 
fore we don't plot them in the figures in order to keep them 
clear. Again, we stress the fact that our definition of resolvable 
source is rather arbitrary, and does not take into account for 
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<5t gw [ns] 

Figure 8. Cumulative number 7V(<5t gw ) of circular binaries (thin 
lines) and bursting triplets (thick lines) emitting over a given 
5tGW threshold as a function of <54qw In each panels the differ- 
ent linestyles refer to the Tu-SA (solid), Tu-DA (long-dashed) and 
Tu-NA (short-dashed) models. The fraction of triplets assumed is 
labeled in each panel. 

the peculiar shape of the burst, we then consider these figures 
as lower limits to the actual detectability of these bursts. 

5.3 Signal samples in the time domain 

To give a feeling of how the actual signals would appear, 
we also computed residuals in the time domain for selected 
sources. To this purpose, we evolved the system using equa- 
tions (27)-(31) of lBarack fc Cutleij (|2004l ) assuming non spin- 
ning SMBHs. We t hen computed all the components h+, n 
and h x ,n (following IPierro. et qZ.lfeoOll ) and we finally evalu- 
ated the residuals R(T) integrating equation (|lip . The actual 
shape of the residuals is rather complex and depends on the 
geometry of the system: the relative orientation of the source 
to the pulsars (encoded in the direction cosines a, j3, and 7 
in equation (|lip ); the polarization angle of the source ^; the 
inclination i; the initial phase of the orbit <E>o; and an angle 
4> p describing t he orientation of the pe riastron in the orbital 
plane (see, e.g., iBarack fc Cutleill2004l for a definition of all 
these quantities). 

Examples of the phenomenology of bursting sources are 
given in figure \W\ for a sample of eccentric systems found in 
one selected realisation of the model Tu-DA. In the left pan- 
els we show the three brightest resolvable sources, while in 
the right panels we show three of the brightest bursts which 
would be unresolvable according to our definition, because 
their power spectra would be overwhelmed by the signal pro- 
duced by the standard circular binaries found in the realiza- 
tion. Parameters of the binaries are given in table [2] Bursts 
can be generated by very eccentric-long period binaries (as in 
the two lower panels), or by relatively short-period systems 
(e.g., central left panel), in which case multiple bursts are vis- 




<5t gw [ns] 

Figure 9. Same as figure [8] but considering only the resolvable 
sources in the computation of N(8t gw ) according to equation (1241 . 
Linestyle as in figure [8] 
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Table 2. Parameters of the sources plotted in figure [T0l Rows in 
the table (from the top to the bottom), correspond to panels of 
figure [T0l considered counterclockwise, starting from the upper left 
panel. 

ible in the observation times. The width of the burst depends 
on the periastron passage timescale: systems with T v <C T b B 
produce narrow features in the data stream (e.g., lower left 
panel), while systems with T p w T b s give a characteristic 
bump shaping all over the data span (e.g. upper right panel). 
Given the integral nature of the signal (equation (JTTJ) ) , its 
shape is also heavily dependent on $0 (e.g. the cumulative 
residual can be positive or negative depending on the binary 
orbital phase at the beginning of the detection), on <f> p and 
on ty. The inclination of the source i and its aperture angle 
to the pulsar, determine the amplitude of the signal. 

5.4 A note for LISA 

We also collected catalogs of systems bursting in the LISA 
window, to check for detectability. Unfortunately, prospects 
for detection with LISA are not as promising as for PTAs. 
In a total of 750 realization of the ten different models, we 
found ~ 50 sources bursting in the LISA window produc- 
ing an SNR> 0.1. Unfortunately, none of them produced 
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Figure 10. Examples of timing residual, found in a particular 
realization of the model Tu-DA, computed according to equation 
WW . Different linestyles correspond to different aperture angle 8 
between the pulsar and the source; 9 = it/2 (dotted), 7r/3 (short- 
dashed), 7r/6 (long-dashed). In all the cases ^ = 7r/4 and i = tt/3. 
<j> v is random and $o is chosen so that the burst occurs during 
the observation. The rms residual computed according to equation 
1 1 41 ) are also shown. Parameters of the sources are listed in table 

m 

a SNR> 8, necessary for a confident detection. We then 
conclude, that even with a consistent population of SMBH 
triplets forming during the cosmic history, burst from mas- 
sive (say, M ~ 10 s Mq) eccentric binaries are unlikely to be 
produced at a significant rate for LISA. On the other hand, 
if formation of triple systems was common in the past, for 
system in the LISA mass range (~ 10 5 — lCFM©), very pe- 
culiar signals from coalescing eccentric binaries may be com- 
mon in the data stream. However, this is beyond the scope 
of the present paper, where we focused on massive binaries 
(M > 1O 7 M ) only. 



6 CONCLUSIONS 

We have addressed in this work three different points in the 
evolution of triplets of SMBHs in the Universe: The Astro- 
dynamics of the system, the potential GW signature and the 
detectability. 

We have performed eight different direct-summation 
N— body simulations, one including more than half a million 
of particles, to calibrate 1,000 3-body scattering experiments, 
which include post-Newtonian corrections, in order to have 
a statistical description of the system. Both numerical tools 
agree that the inner binary of SMBHs will go through a phase 
of extremely high eccentricity, which is the motivation for the 
rest of the work. 

These three-body excitations of episodic high eccentric- 
ity configurations of the close SMBH binary produce inter- 
esting GW bursts that may be detectable with forthcoming 



experiments such as PTAs and LISA. The extreme eccentrici- 
ties of such bursts on one hand would leave a very distinctive 
signature, but on the other require the development of ap- 
propriate analysis techniques. 

To compute likely event rates, we extracted catalogues 
of merging galaxies from the Millennium Run, and we pop- 
ulated them with SMBHs following the known MBH-bulges 
relations. We then estimated the fractions of triplets and their 
eccentricity distribution and we computed the induced signals 
in both PTAs and the LISA detector. 

We found that, depending on the details of the SMBH 
population model, if the fraction of triplets is ^ 0.1, few to a 
hundred of GW bursts would be produced at a > 1 ns level 
in the PTA frequency domain. Most of the signals will be 
washed out in the confusion noise due to the emission of 'ordi- 
nary' low eccentric binaries. However, their peculiar features 
may guide the development of targeted data analysis tech- 
niques, that may help to recognize them even if overwhelmed 
by the confusion noise. Employing a minimal criterion for 
source resolvability (which provides a strict lower limit), we 
found that less than one system may be actually pinned down 
at ns precisions. By running several dozens of Monte Carlo 
realization of the signal from the cosmological population of 
SMBH binaries and triplets we quantified a statistical 2—50% 
chance of having a resolvable burst in the Universe (assum- 
ing 10 yrs of observation). The probability for detection with 
LISA is essentially nil. However, we stress the fact that we 
focused on systems with M > 10 7 Mq; our results then sim- 
ply imply that it is extremely unlikely that a system which 
would normally emit outside the LISA range will produce a 
burst in the LISA window because of resonant three body in- 
teractions. On the other hand, if a consistent fraction of light 
binaries (M < 10 7 Mq) is involved in triple systems, we may 
expect several eccentricity-driven coalescences to be observed 
by LISA. This eventuality would call for the development of 
extremely eccentric templates (e > 0.9) for merging SMBHs, 
and of adequate analysis techniques to extract the signal. 
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